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. At high energy, cosmic rays can only be studied by measuring the extensive air showers 

' they produce in the atmosphere of the Earth. Although the main features of air showers can 

(N ■ be understood within a simple model of successive interactions, detailed simulations and 

■ a realistic description of particle production are needed to calculate observables relevant 
^ ' to air shower experiments. Currently hadronic interaction models are the main source of 

, uncertainty of such simulations. We will study the effect of using different hadronic models 

QQ ■ available in CORSIKA and CONEX on extensive air shower predictions. 

^ I PACS: 96.50.sb 96.50.sd 

^ , Key words: cosmic rays, air showers, hadronic interaction models 

o : 
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, 1 Introduction 

(N ■ 

' Due to the steeply falling energy spectrum of cosmic rays, direct detection by 

■ satellite- or balloon-borne instruments is only possible up to about ~ 10^^ eV. 
I Fortunately, at such high energy, the cascades of secondary particles produced by 
. cosmic rays reach the ground and can be detected in coincidence experiments. The 

Oh' cascades are called extensive air showers and are routinely used to make indirect 

Q . measurements of high energy cosmic rays. The upper limit of the detectable energy 

i—^ ' is given by the area and exposure time of the detector. For instance, the Pierre 

c/3 . Auger Observatory J] , whose Southern detector is currently under construction in 

' Argentina, is designed to detect particles of ^ 10^° eV for which the flux is less 

^ . than one particle per km^ and century. 

Air showers can be observed using different detection techniques. The most 
frequently employed technique is the measurement of secondary particles reaching 
' ground. Using an array of particle detectors (for example, sensitive to and n^), 

the arrival direction and information on mass and energy of the primary cosmic 
ray can be reconstructed. The main observables are the number and the lateral and 
temporal distributions of the different secondary particles. At energies above ^ 10^^ 
eV, the longitudinal profile of a shower can be directly observed by measuring the 
fluorescence light induced by the charged particles traversing the atmosphere. Two 
main observables can be extracted from the longitudinal shower profile: the energy 
deposit or the number of particles, -/Vmax, at the shower maximum and Xmax, the 
atmospheric depth of the maximum. Again, these quantities can be used to estimate 
the energy and mass of the primary particles. Shower-to-shower fluctuations of all 
observables provide also very useful composition information. 

As a consequence of the indirect character of the measurement, detailed sim- 
ulations of air showers are needed to extract information on the primary particle 
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from shower observables. Whereas electromagnetic interactions are well understood 
within perturbative QED, hadronic multiparticle production cannot be calculated 
within QCD from first principles. Differences in modelling hadronic interactions, 
which cannot be resolved by current accelerator data, are the main source of un- 
certainty of air shower predictions |21IS1- 

In this article, we will discuss the relation between hadronic multiparticle pro- 
duction and air shower observables. In Sec. [3 a very simple toy model, based on 
the model by Heitler 4 from the fifties, is used to introduce the basic features of 
air showers and their relation to multiparticle production. Using detailed Monte 
Carlo simulations done with CORSIKA |S] and CONEX 6 the uncertainty of cosmic 
ray measurements due to hadronic interaction models is illustrated in Sec. |31 Here 
we concentrate on the uncertainties implied by different models. In addition, each 
model has a number of parameters that influence the extrapolation to high energy. 
Studies of the uncertainty range due to variation of model parameters within a 
single model can be found in |7| |S| El 021 • In general, the differences between mod- 
els are comparable or bigger than the differences of the predictions that can be 
obtained within a single model by parameter adjustment. 



2 Heitler's Model 



To qualitatively describe the dependence of shower development on some basic 
parameters of particle interaction, decay and production, a very simple toy model 
can be used. Although initially developed for electromagnetic (EM) showers it 
can also be applied to hadronic showers [TTl. 



2.1 Electromagnetic showers 

For simplicity, instead of having three particle types (7, e"*" and e~) like in elec- 
tromagnetic (EM) showers, we will consider only one particle with energy E with 
only one EM interaction producing two new particles with energy E/2 after a fixed 
interaction length of Ae, see Fig.^a). 

Denoting with n the number of generations (consecutive interactions), the num- 
ber of particles at a given depth X — n ■ follows from 

iV(X) = 2" = 2-^/^% (1) 

with the energy E per particle for a given primary energy Eq being 

^(^) = (2) 

Defining the critical energy Ec (~ 85 MeV in air) as the energy below which 
energy loss processes dominate over particle production, one can make the assump- 
tion that the shower maximum is reached at a depth at which the energy of the 
secondary particles reaches Ef.- Then two main shower observables are given by 

^max = ^ and X^^,^iEo)^Xe-\n(^^y (3) 
2 A Czech. J. Phys. 56 (2005) 
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Fig. 1. Schematic view of (a) electromagnetic and (b) hadronic cascades. In the latter, 
dashed lines represent neutral particles (vr") and solid lines charged particles (tt*). Only 
one charged hadron interaction is shown for each generation. 

This simplified picture does not reproduce the detailed behavior of an EM shower, 
but two important features are well described: the number of particles at shower 
maximum is proportional to Eq and the depth of shower maximum depends loga- 
rithmically on the primary energy Eq. 

2.2 Hadronic showers 

Generalizing this idea, a hadronic interaction of a particle with energy E is assumed 
to produce ntot new particles with energy E/ntot, two third of which being charged 
particles rich (charged pions) and one third being neutral particles nnout (neutral 
pions), as shown Fig.^b). Neutral particles decay immediately into EM particles 
(tt^ —^ 27) . After having travelled a distance corresponding to the mean interaction 
length Aino, charged particles re-interact with air nuclei if their energy is greater 
than some typical decay energy i^dec- 

Energy transfer 

In each hadronic interaction, one third of the energy is transferred to the EM 
shower component. After n generations the energy in the hadronic and EM com- 
ponents is given by 



Ehad = 
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Eem 



Eo, (5) 



where n will be calculated later. Simulations show that the number of generations 
is typically about 5 to 6 jl2| . 

Even in an air shower initiated by a hadron, most of the energy is carried by EM 
particles 90% for n = 6). Hence the depth of shower maximum is given by that 
of the EM shower component, X^^^. As the first hadronic interaction produces EM 
particles of energy ~ Eg /ntot one gets 

XmaxiEo) ^ Mnc + X^^^^{Eo/ntot) (6) 

where X-mc is the hadronic interaction length. This simplified expression for the 
depth of maximum neglects the EM sub-showers initiated by hadrons of later gen- 
erations. The inclusion of higher hadronic generations does not change the structure 
of Eq. 0, only the coefficients change (see, for example, 13 ). 

Muon component 

To keep the picture simple, we assume that all charged hadrons decay into 
muons when their energy reaches -Edec- In a real shower, this limit can be seen 
as the characteristic energy where interaction length and decay length of charged 
pions are similar (about 150 GeV for pions). By construction, charged particles will 
reach the energy i^dcc after n interactions 

^--(^- 

Since one muon is produced in the decay of each charged particle, we get 

^.-<h-fi^y, (9) 

with a — In rich/ hmtot ~ 0-82 . . . 0.95 |13| . The number of muons produced in an 
air shower depends not only on the primary energy and air density, but also on the 
charged and total particle multiplicities of hadronic interactions. 

It should be kept in mind that the parameters of the model are only effective 
quantities and are not identical to the respective quantities measured at accel- 
erators. In particular, the approximation of all secondary particles carrying the 
same energy is only motivated by the fact that it allows us to obtain simple, closed 
expressions. The well-known leading particle effect, typically quantified by the (in)- 
elasticity of an interaction, can be implemented in the model [IV but will not be 
considered here. 
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Superposition model 

In case of a nucleus being the primary particle, one can use the superposi- 
tion model to deduce the main observables from the above-written formulas. In 
this model, a nucleus with mass A and energy Eq is considered as A independent 
nucleons with energy Ei^ = Eq/A. This leads to 



Note that, whereas there is no mass dependence of the number of charged particles 
at shower maximum, both the number of muons and the depth of maximum depend 
on the mass of the primary particle. The heavier the shower-initiating particle the 
more muons are expected for a given primary energy. 



It is clear that the model described above is only giving a very much over- 
simplified account of air shower physics. However, the model allows us to qualita- 
tively understand the dependence of many air shower observables on the character- 
istics of hadronic particle production. The parameters of hadron production being 
most important for air shower development are the cross section (or mean free path 
Aine), the multiplicity of secondary particles of high energy, ntot, and the produc- 
tion ratio of neutral to charged particles. Unfortunately, these parameters are not 
well constrained by particle production measurements at accelerators. Depending 
on the assumptions on how to extrapolate existing accelerator data, the predictions 
of hadronic interaction models differ considerably. 

3.1 Hadronic interaction models 

There are several hadronic interaction models commonly used to simulate air show- 
ers. For high energy interactions (-Eiab ^ 100 GeV), these models are DPMJET 11.55 
and III [HI, NeXus 2 and 3.97 ,15, 16j, QGSJET 01 and II 17, 18, 19 , and SIBYLL 2.1 
ISnilSES- The physics models and assumptions are discussed in, for example, pTij . 
All the high-energy interaction models reproduce accelerator data reasonably well 
but predict different extrapolations above E'cms ~1.8 TeV (i?iab ~ 10^^ eV) that 
lead to very different results at high energy |Sl|2ni- The situation is different at 
low energy where several measurements from fixed target experiments are avail- 
able |12j . There one of the main problems is the extrapolation of measurements to 
the very forward phase space region close to the beam direction and the lack of 
measurements of pion-induced interactions. At low energy, models based on data 
parameterization and/or microscopic models such as fluka [23 j GHEISHA j25,, or 
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UrQMD 1.3 are used. Differences in the results of those models mainly play a 
role for the number of muons far from the shower axis [271 128j . 

In the following we will only focus on high energy interaction models. A first 
illustration of differences of air shower predictions due to model uncertainties is 
shown in Fig. [21 The shower maximum is shifted from one model to another by up 
to 40 g/cm^. The difference is bigger for proton induced showers than that of iron. 
The small differences of the predictions of iron showers can partially be explained 
by the superposition model. Iron showers can be considered as many proton showers 
of 1/56 of the primary energy since. At lower energy, closer to existing collider data, 
the model predictions are more similar. On the other hand the number of particles 
at shower maximum is almost model-independent. 




Fig. 2. Mean longitudinal profile of vertical iron, proton and 7 induced showers at 10 eV 
using diflferent high-energy hadronic interaction models (100 showers averaged). 

3.2 Simulation results 

Using the air shower simulation packages CORSIKA and CONEX and various high- 
energy hadronic interaction models, we can test the model dependence of different 
energy and mass estimators of air shower experiments and get an estimate of the 
resulting uncertainties. 

Energy estimation 

In case of a fluorescence light detector, the energy is measured by integrating 
the signal of the observed longitudinal ionization energy deposit profile. To account 
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for energy carried by neutrinos and partially also muons, one has to convert the 
observed, calorimetric energy to the total shower energy. The conversion factor is 
energy- and mass-dependent. Predictions for the conversion factor obtained from 
simulations done with CONEX are shown in Fig. |3| From Eq. Q , we know that 



1.3 



1.25 



1.2 



1.15 



1.1 



1.05 



Fe 



• QGSJET 01 
■ SIBYLL 2.1 
A NEXUS 3.97 

— mean Mass 

- - - mean Proton 
mean Iron 




17 



17.5 



18 



18.5 



19 



19.5 



20 



20.5 



log (E /eV) 

*'10 tot 



Fig. 3. Mean factor for the conversion of observed (calorimetric) energy to total energy for 
iron (dotted) or proton (dashed) induced showers |29| . The conversion factor is shown for 
QGSJET 01 (circles), SIBYLL 2.1 (squares) or neXus 3.97 (triangles). The mean conversion 
factor is calculated by averaging all proton and iron predictions. 

this conversion factor depends mainly on the effective number of generations, which 
in turn is only logarithmically dependent on ntot- Therefore, for a given mass, the 
uncertainty due to the hadronic model is only about 2%. Furthermore, the unknown 
mass of the primary particles leads to an uncertainty of less than 3% on the energy 
estimation above 10^^ eV. 

Ground based detectors use the density of particles at a given distance from 
the shower core to reconstruct the shower energy. In case of the Pierre Auger 
Observatory, the Cherenkov density in units of that of vertical muons (vertical 
equivalent muon, 1 VEM « 240 MeV energy deposit) is measured at a lateral 
distance of 1000 m. CORSIKA simulations of this quantity are shown for vertical 
proton and iron induced showers at lO^^eV in Fig. 21 In this case, the uncertainty 
between models for proton induced showers is about ±7%. If one takes into account 
the mass uncertainty, the theoretical uncertainty increases to about ±17%. At least 
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qualitatively, the difference between QGSJET and SIBYLL can be understood in terms 
of the a parameter of the Heitler-type model. The a parameter is lower for the 
SIBYLL model than for QGSJET, 0.84 and 0.85, respectively fTT. This small difference 
is sufficient to give a 15% difference in the number of muons which contribute about 
50% to the Cherenkov density. 
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Fig. 4. Mean lateral distribution function of Cherenkov density for 10 eV vertical proton 
and iron induced showers and different high-energy hadronic interaction models. 



Mass estimation 

It is desirable to estimate both the energy and the mass of the primary particles 
at the same time. First of all the composition of cosmic rays holds important clues 
on the nature of their sources propagation processes. Furthermore one can improve 
the energy reconstruction considerably by knowing the mass of the primary particle. 

In an experiment employing the fluorescence light technique, the mass of the 
primary particle can be estimated through the measurement of the depth of shower 
maximum. In Fig. [3 CONEX simulations of the mean depth of shower maximum, 
(^max), are shown together with a compilation of measurements. The model depen- 
dence of (Xiiiax) is obvious. The mean depth of shower maximum mainly depends 
on the hadronic inelastic cross section, the multiplicity as well as the inelasticity, 
see Eq. (|7||. For example, neXus (dash-dotted line with stars) is characterized by 
the highest inelastic cross section and the lowest multiplicity. Therefore it is not 
surprising that it predicts the highest (Xmax)- On the other hand, QGSJET 01 (full 
line with open circles) has the lowest inelastic cross section and the highest multi- 
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Fig. 5. Mean X^ax for proton and iron induced showers as a function of the primary 
energy. Predictions of different high-energy hadronic interaction models are compared to 
data. Refs. to the data can be found in 1301. 



plicity, which explains why its (Xmax) is small. The difference between the model 
predictions is limiting the conclusions that can be drawn from the measurements. 
For example, in a simple interpretation of the experimental data at the highest 
energy, using QGSJET 01 will lead to a composition dominated by protons rather 
than a mixed composition as would be obtained if neXus is taken as reference. 

The mass dependence of the correlation of the charged particle and muon mul- 
tiplicities of showers is used in surface array experiments for both energy and com- 
position reconstruction. According to the prediction of the superposition model, see 
Eqs. (|1U|I and (|12|l . one would expect a separation of proton and iron induced show- 
ers at fixed energy in the Nch — plane. However, the sensitivity of this method 
of mass reconstruction is limited by its model dependence and shower-to-shower 
fluctuations. In Fig.El we show the charged particle and muon shower sizes of 50 
proton and iron showers at Auger observation level (1450m a.s.l.). Each point corre- 
sponds to a 30° inclined shower simulated at 10^° eV. For clarity, the predictions of 
only two models are shown, QGSJET 01 (triangles) and neXus (squares). There is 
a clear separation between proton (full symbols) and iron (open symbols) initiated 
showers for each model. The number of electrons at ground is slightly higher for 
proton induced showers because A^max is larger and the maximum is closer to the 
ground. As expected, the number of muons in iron induced showers is higher than 
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in proton ones by about a factor of 1.4 which corresponds to A^^" with a ^ 0.92. 
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Fig. 6. Correlation of charged particle and muon multiplicites of proton (filled symbols) 
and iron (open symbols) induced showers at 10^" eV. Each point represents an individual 
shower simulated with CONEX. Triangles show results obtained with QGSJET 01 and squares 

that of neXus 3.97. 



Applying another model, for instance NeXus 3.97 (sibyll gives similar results in 
this case), we notice that proton and iron initiated showers are even more separated 
than in QGSJET simulations (a ^ 0.86 is lower and the difference between (Xmax) 
for proton and iron is larger). However, iron induced showers from NeXus populate 
a region in the plot that overlaps that of proton induced showers from QGSJET 01. 
This model dependence makes it very difficult to interpret experimental data. Using 
not only the mean multiplicities but also the shower-to-shower fluctuations and 
correlations in the analysis, one can try to estimate the mass within a given hadronic 
interaction model ,31. , and even exclude models if the data can not be reproduced at 
all. Finally it should be noted that the fluctuations and also the model dependence 
are somewhat reduced if particle densities at about 600 to 1000m from the shower 
core are considered instead of total multiplicity (see, for example, IS^ITU] ). 

Finally it should be mentioned that there are a number of alternative high 
energy extrapolations of hadronic multiparticle production suggested in literature 
(for example, see [331|^|3S]) that have not been discussed in this work. It can be 
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expected that the uncertainty of the composition reconstruction increases consider- 
ably if these scenarios are included. On the other hand, the energy reconstruction 
by means of fluorescence light detection is very robust and will change only slightly. 



4 Summary 

Using a simple cascade model, it is possible to find the main parameters of 
hadronic interactions that influence air shower predictions. These parameters, name- 
ly the inelastic cross sections, the secondary particle multiplicity, the inelasticity, 
and the ratio of charged to neutral hadrons, depend of the hadronic interaction 
model. As a consequence, realistic simulations of hadron induced air-showers are 
model-dependent, leading to theoretical uncertainties in the analysis of experimen- 
tal data. For a ground based detector, the model-related systematic error on energy 
estimation can be as large as 17% at 10^^ eV if the mass of the primary particle 
is unknown. The theoretical uncertainties of the energy reconstruction are much 
smaller for fluorescence light detectors (less than 5% even for unknown primary 
particle mass). The model dependence of the primary mass estimation is crucial 
and currently the mass composition can only be derived for a given hadronic model. 
On the other hand, cosmic ray data of multi-observable detectors can be used to 
test hadronic interaction models at energies much higher than those reached in 
laboratory. 
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